/*******************************************************************************
 2HDMC - two-Higgs-doublet model calculator
 Demo program

 http://2hdmc.hepforge.org
*******************************************************************************/
#include "Constraints.h"
#include "DecayTable.h"
#include "HBHS.h"
#include "SM.h"
#include "THDM.h"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>

using namespace std;

bool isFileExists(string &name){
	ifstream f(name.c_str());
	return f.good();
}

int main(int argc, char *argv[]) {

	if(argc!=12){

		std::cout<< "Please use it as: ./single_point mh mH mC mA sba tb l6 l7 m12_2 2HDM_type output_file_name" <<std::endl;
		exit(0);

	}
	string file_name = argv[11];
	//fstream f;
	FILE *fp; 
	if (isFileExists(file_name)){
		fp = fopen(file_name.c_str(),"a+");
	}
	else{
		fp = fopen(file_name.c_str(),"w+");
		fprintf(fp, "%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s\
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s %12s %12s %12s %12s %12s %12s \
%12s %12s %12s \
%12s %12s %12s %12s \n",
				"mh", "mH", "mA", "mC" , "sba", "tb", "l6", "l7", "m12_2", "type",  
				"stability", "unitarity", "perturb",  "hb_result[0]", "hb_result[1]", "hb_result[2]", "hb_result[3]", 
				"hb_chan[0]", "hb_chan[1]", "hb_chan[2]", "hb_chan[3]", "hb_obsr[0]", "hb_obsr[1]", "hb_obsr[2]", "hb_obsr[3]", "hs_chisq", "hs_pvalue",
				"S", "T", "U", "V", "W", "X", "Delta_rho", "Delta_amu", 
				"Tot_Ga(h)", "Ga(hWW)","Ga(hZZ)","Ga(hmumu)","Ga(htata)", "Ga(hss)", "Ga(hbb)", "Ga(huu)", "Ga(hcc)", "Ga(hgg)", "Ga(hgaga)", "Ga(hgaZ)", "Ga(hAA)", "Ga(hHpHm)","Ga(hhh)", "Ga(hZA)", "Ga(hWHp)",  
				             "Br(hWW)","Br(hZZ)","Br(hmumu)","Br(htata)", "Br(hss)", "Br(hbb)", "Br(huu)", "Br(hcc)", "Br(hgg)", "Br(hgaga)", "Br(hgaZ)", "Br(hAA)", "Br(hHpHm)","Br(hhh)", "Br(hZA)", "Br(hWHp)",  
				"Tot_Ga(H)", "Ga(HWW)","Ga(HZZ)","Ga(Hmumu)","Ga(Htata)", "Ga(Hss)", "Ga(Hbb)", "Ga(Huu)", "Ga(Hcc)", "Ga(Hgg)", "Ga(Hgaga)", "Ga(HgaZ)", "Ga(HAA)", "Ga(HHpHm)","Ga(Hhh)", "Ga(HZA)", "Ga(HWHp)",  
				             "Br(HWW)","Br(HZZ)","Br(Hmumu)","Br(Htata)", "Br(Hss)", "Br(Hbb)", "Br(Huu)", "Br(Hcc)", "Br(Hgg)", "Br(Hgaga)", "Br(HgaZ)", "Br(HAA)", "Br(HHpHm)","Br(Hhh)", "Br(HZA)", "Br(HWHp)",  
				"Tot_Ga(A)", "Ga(AWW)","Ga(AZZ)","Ga(Amumu)","Ga(Atata)", "Ga(Ass)", "Ga(Abb)", "Ga(Auu)", "Ga(Acc)", "Ga(Agg)", "Ga(Agaga)", "Ga(AgaZ)", "Ga(AZh)", "Ga(AZH)", "Ga(AWHp)",  
				             "Br(AWW)","Br(AZZ)","Br(Amumu)","Br(Atata)", "Br(Ass)", "Br(Abb)", "Br(Auu)", "Br(Acc)", "Br(Agg)", "Br(Agaga)", "Br(AgaZ)", "Br(AZh)", "Br(AZH)", "Br(AWHp)",  
				"Tot_Ga(Hp)","Ga(Hptb)","Ga(Hpcs)","Ga(Hpcb)","Ga(Hptanu)", "Ga(Hpmunu)", "Ga(Hpenu)", "Ga(HpWh)", "Ga(HpWH)", "Ga(HpWA)", 
				             "Br(Hptb)","Br(Hpcs)","Br(Hpcb)","Br(Hptanu)", "Br(Hpmunu)", "Br(Hpenu)", "Br(HpWh)", "Br(HpWH)", "Br(HpWA)", 
				"Tot_Ga(top)","Ga(tHpb)","Br(tHpb)",
				"CS_signalH", "CS_Hhh", "CS_bbgaga", "CS_tHpb");
	}
	// Reference SM Higgs mass for EW precision observables
	// Create SM and set parameters
	SM sm;
	sm.set_qmass_pole(6, 172.5);
	sm.set_qmass_pole(5, 4.75);
	sm.set_qmass_pole(4, 1.42);
	sm.set_lmass_pole(3, 1.77684);
	sm.set_alpha(1. / 127.934);
	sm.set_alpha0(1. / 137.0359997);
	sm.set_alpha_s(0.119);
	sm.set_MZ(91.15349);
	sm.set_MW(80.36951);
	sm.set_gamma_Z(2.49581);
	sm.set_gamma_W(2.08856);
	sm.set_GF(1.16637E-5);
	
	// Create 2HDM and set SM parameters
	THDM model;
	model.set_SM(sm);
	
	// Set parameters of the 2HDM in the 'physical' basis
	double mh       = atof(argv[1]);
	double mH       = atof(argv[2]);
	double mA       = atof(argv[3]);
	double mC       = atof(argv[4]);
	double sba      = atof(argv[5]);
	double tb       = atof(argv[6]);
	double lambda_6 = atof(argv[7]);
	double lambda_7 = atof(argv[8]);
	double m12_2    = atof(argv[9]);
	int    typ      = atoi(argv[10]);
	
	//std::cout << typ <<std::endl;
	bool pset = model.set_param_phys(mh, mH, mA, mC, sba, lambda_6, lambda_7, m12_2, tb);
	
	if (!pset) {
		cerr << "The specified parameters are not valid\n";
		return -1;
	}

	// Set Yukawa couplings to 2HDM-type 
	model.set_yukawas_type(typ);
	
	// Prepare to calculate decay widths
	DecayTable table(model);
	
  //############  h decays #####################################
  double gamma_total_h = table.get_gammatot_h(1);
  double gamma_hWW     = table.get_gamma_hvv(1,3)  ;
  double gamma_hZZ     = table.get_gamma_hvv(1,2)  ;
  double gamma_hmumu   = table.get_gamma_hll(1,2,2);
  double gamma_htata   = table.get_gamma_hll(1,3,3);
  double gamma_hss     = table.get_gamma_hdd(1,2,2);
  double gamma_hbb     = table.get_gamma_hdd(1,3,3);
  double gamma_huu     = table.get_gamma_huu(1,1,1);
  double gamma_hcc     = table.get_gamma_huu(1,2,2);
  double gamma_hgg     = table.get_gamma_hgg(1)    ;
  double gamma_hgaga   = table.get_gamma_hgaga(1)  ;
  double gamma_hgaZ    = table.get_gamma_hZga(1)   ;
  double gamma_hAA     = table.get_gamma_hhh(1,3,3);
  double gamma_hHpHm   = table.get_gamma_hhh(1,4,4);
  double gamma_hhh     = table.get_gamma_hhh(1,1,1);
  double gamma_hZA     = table.get_gamma_hvh(1,2,3);
  double gamma_hWHp    = table.get_gamma_hvh(1,4,4);

  double Br_hWW     = gamma_hWW   / gamma_total_h;
  double Br_hZZ     = gamma_hZZ   / gamma_total_h;
  double Br_hmumu   = gamma_hmumu / gamma_total_h;
  double Br_htata   = gamma_htata / gamma_total_h;
  double Br_hss     = gamma_hss   / gamma_total_h;
  double Br_hbb     = gamma_hbb   / gamma_total_h;
  double Br_huu     = gamma_huu   / gamma_total_h;
  double Br_hcc     = gamma_hcc   / gamma_total_h;
  double Br_hgg     = gamma_hgg   / gamma_total_h;
  double Br_hgaga   = gamma_hgaga / gamma_total_h;
  double Br_hgaZ    = gamma_hgaZ  / gamma_total_h;
  double Br_hAA     = gamma_hAA   / gamma_total_h;
  double Br_hHpHm   = gamma_hHpHm / gamma_total_h;
  double Br_hhh     = gamma_hhh   / gamma_total_h;
  double Br_hZA     = gamma_hZA   / gamma_total_h;
  double Br_hWHp    = gamma_hWHp  / gamma_total_h;

  // H decays #####################################################
  double gamma_total_H = table.get_gammatot_h(2);
  double gamma_HWW     = table.get_gamma_hvv(2,3)  ;
  double gamma_HZZ     = table.get_gamma_hvv(2,2)  ;
  double gamma_Hmumu   = table.get_gamma_hll(2,2,2);
  double gamma_Htata   = table.get_gamma_hll(2,3,3);
  double gamma_Hss     = table.get_gamma_hdd(2,2,2);
  double gamma_Hbb     = table.get_gamma_hdd(2,3,3);
  double gamma_Huu     = table.get_gamma_huu(2,1,1);
  double gamma_Hcc     = table.get_gamma_huu(2,2,2);
  double gamma_Hgg     = table.get_gamma_hgg(2)    ;
  double gamma_Hgaga   = table.get_gamma_hgaga(2)  ;
  double gamma_HgaZ    = table.get_gamma_hZga(2)   ;
  double gamma_HAA     = table.get_gamma_hhh(2,3,3);
  double gamma_HHpHm   = table.get_gamma_hhh(2,4,4);
  double gamma_Hhh     = table.get_gamma_hhh(2,1,1);
  double gamma_HZA     = table.get_gamma_hvh(2,2,3);
  double gamma_HWHp    = table.get_gamma_hvh(2,4,4);

  double Br_HWW     = gamma_HWW   / gamma_total_H;
  double Br_HZZ     = gamma_HZZ   / gamma_total_H;
  double Br_Hmumu   = gamma_Hmumu / gamma_total_H;
  double Br_Htata   = gamma_Htata / gamma_total_H;
  double Br_Hss     = gamma_Hss   / gamma_total_H;
  double Br_Hbb     = gamma_Hbb   / gamma_total_H;
  double Br_Huu     = gamma_Huu   / gamma_total_H;
  double Br_Hcc     = gamma_Hcc   / gamma_total_H;
  double Br_Hgg     = gamma_Hgg   / gamma_total_H;
  double Br_Hgaga   = gamma_Hgaga / gamma_total_H;
  double Br_HgaZ    = gamma_HgaZ  / gamma_total_H;
  double Br_HAA     = gamma_HAA   / gamma_total_H;
  double Br_HHpHm   = gamma_HHpHm / gamma_total_H;
  double Br_Hhh     = gamma_Hhh   / gamma_total_H;
  double Br_HZA     = gamma_HZA   / gamma_total_H;
  double Br_HWHp    = gamma_HWHp  / gamma_total_H;
  
  // A decays #####################################################
  double gamma_total_A = table.get_gammatot_h(3);
  double gamma_AWW     = table.get_gamma_hvv(3,3)  ;
  double gamma_AZZ     = table.get_gamma_hvv(3,2)  ;
  double gamma_Amumu   = table.get_gamma_hll(3,2,2);
  double gamma_Atata   = table.get_gamma_hll(3,3,3);
  double gamma_Ass     = table.get_gamma_hdd(3,2,2);
  double gamma_Abb     = table.get_gamma_hdd(3,3,3);
  double gamma_Auu     = table.get_gamma_huu(3,1,1);
  double gamma_Acc     = table.get_gamma_huu(3,2,2);
  double gamma_Agg     = table.get_gamma_hgg(3)    ;
  double gamma_Agaga   = table.get_gamma_hgaga(3)  ;
  double gamma_AgaZ    = table.get_gamma_hZga(3)   ;
  double gamma_AZh     = table.get_gamma_hvh(3,2,1);
  double gamma_AZH     = table.get_gamma_hvh(3,2,2);
  double gamma_AWHp    = table.get_gamma_hvh(3,4,4);

  double Br_AWW     = gamma_AWW   / gamma_total_A;
  double Br_AZZ     = gamma_AZZ   / gamma_total_A;
  double Br_Amumu   = gamma_Amumu / gamma_total_A;
  double Br_Atata   = gamma_Atata / gamma_total_A;
  double Br_Ass     = gamma_Ass   / gamma_total_A;
  double Br_Abb     = gamma_Abb   / gamma_total_A;
  double Br_Auu     = gamma_Auu   / gamma_total_A;
  double Br_Acc     = gamma_Acc   / gamma_total_A;
  double Br_Agg     = gamma_Agg   / gamma_total_A;
  double Br_Agaga   = gamma_Agaga / gamma_total_A;
  double Br_AgaZ    = gamma_AgaZ  / gamma_total_A;
  double Br_AZh     = gamma_AZh   / gamma_total_A;
  double Br_AZH     = gamma_AZH   / gamma_total_A;
  double Br_AWHp    = gamma_AWHp  / gamma_total_A;
  
 
 // hp decays #####################################################  
  double gamma_total_Hp = table.get_gammatot_h(4);
  double gamma_Hptb   = table.get_gamma_hdu(4,3,3);
  double gamma_Hpcs   = table.get_gamma_hdu(4,2,2);
  double gamma_Hpcb   = table.get_gamma_hdu(4,3,2); 
  double gamma_Hptanu = table.get_gamma_hln(4,3,3); 
  double gamma_Hpmunu = table.get_gamma_hln(4,2,2);   
  double gamma_Hpenu  = table.get_gamma_hln(4,1,1);
  double gamma_HpWh   = table.get_gamma_hvh(4,3,1);
  double gamma_HpWH   = table.get_gamma_hvh(4,3,2);
  double gamma_HpWA   = table.get_gamma_hvh(4,3,3);

  double Br_Hptb   = gamma_Hptb   / gamma_total_Hp;
  double Br_Hpcs   = gamma_Hpcs   / gamma_total_Hp;
  double Br_Hpcb   = gamma_Hpcb   / gamma_total_Hp; 
  double Br_Hptanu = gamma_Hptanu / gamma_total_Hp; 
  double Br_Hpmunu = gamma_Hpmunu / gamma_total_Hp;   
  double Br_Hpenu  = gamma_Hpenu  / gamma_total_Hp;
  double Br_HpWh   = gamma_HpWh   / gamma_total_Hp;
  double Br_HpWH   = gamma_HpWH   / gamma_total_Hp;
  double Br_HpWA   = gamma_HpWA   / gamma_total_Hp;

 // top decays #####################################################
  double gamma_total_top = table.get_gammatot_top();
  double gamma_tHpb   = table.get_gamma_uhd(3,4,3) ; //top -> Hp + b
  double Br_tHpb   = gamma_tHpb / gamma_total_top;

	//const HBHSResult *hbhsres_ptr = nullptr;
	HBHS hbhs{};
	hbhs.HBHS_setup_output_level(0);

	
	//Prepare to calculate observables
	Constraints constr(model);
	
	double S, T, U, V, W, X, mh_ref;
	
	//mh/mH, which is closed to 125.0, is set as SM Higgs
	if(abs(mh-125.0)>abs(mH-125.0)){
		mh_ref = mH;
	}
	else{
		mh_ref = mh;
	}
	constr.oblique_param(mh_ref, S, T, U, V, W, X);
    double rho  = constr.delta_rho(mh_ref);
    double damu = constr.delta_amu()*pow(10,10); //inlarge for 10^10
	
	bool stability= constr.check_stability();
	bool unitarity= constr.check_unitarity();
	bool perturbativity= constr.check_perturbativity();

	const HBHSResult hbhs_result = hbhs.check(model);

    double CS_singleH = hbhs_result.hb.CS_LHC13[1].find("singleH")->second;
	double CS_Hhh     = CS_singleH * Br_Hhh;
	double CS_bbgaga  = CS_Hhh * Br_hbb * Br_hgaga;
    double CS_tHpb = 966.01*Br_tHpb;

	fprintf(fp, "%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12i \
%12i %12i %12i %12i %12i %12i %12i \
%12i %12i %12i %12i %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e\
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e \
%12.5e %12.5e %12.5e %12.5e \n",  
		         mh, mH, mA, mC, sba, tb, lambda_6, lambda_7, m12_2, typ, 
		         stability, unitarity, perturbativity, hbhs_result.hb.result[0], hbhs_result.hb.result[1], hbhs_result.hb.result[2], hbhs_result.hb.result[3], 
		         hbhs_result.hb.channel[0], hbhs_result.hb.channel[1], hbhs_result.hb.channel[2], hbhs_result.hb.channel[3], hbhs_result.hb.obsratio[0], hbhs_result.hb.obsratio[1], hbhs_result.hb.obsratio[2], hbhs_result.hb.obsratio[3], hbhs_result.hs.chisq, hbhs_result.hs.pvalue, 
		         S, T, U, V, W, X, rho, damu,
				 gamma_total_h, gamma_hWW, gamma_hZZ, gamma_hmumu, gamma_htata, gamma_hss, gamma_hbb, gamma_huu, gamma_hcc, gamma_hgg, gamma_hgaga, gamma_hgaZ, gamma_hAA, gamma_hHpHm, gamma_hhh, gamma_hZA, gamma_hWHp,
				                Br_hWW, Br_hZZ, Br_hmumu, Br_htata, Br_hss, Br_hbb, Br_huu, Br_hcc, Br_hgg, Br_hgaga, Br_hgaZ, Br_hAA, Br_hHpHm, Br_hhh, Br_hZA, Br_hWHp,
				 gamma_total_H, gamma_HWW, gamma_HZZ, gamma_Hmumu, gamma_Htata, gamma_Hss, gamma_Hbb, gamma_Huu, gamma_Hcc, gamma_Hgg, gamma_Hgaga, gamma_HgaZ, gamma_HAA, gamma_HHpHm, gamma_Hhh, gamma_HZA, gamma_HWHp,
				                Br_HWW, Br_HZZ, Br_Hmumu, Br_Htata, Br_Hss, Br_Hbb, Br_Huu, Br_Hcc, Br_Hgg, Br_Hgaga, Br_HgaZ, Br_HAA, Br_HHpHm, Br_Hhh, Br_HZA, Br_HWHp,
				 gamma_total_A, gamma_AWW, gamma_AZZ, gamma_Amumu, gamma_Atata, gamma_Ass, gamma_Abb, gamma_Auu, gamma_Acc, gamma_Agg, gamma_Agaga, gamma_AgaZ, gamma_AZh, gamma_AZH, gamma_AWHp,
				                Br_AWW, Br_AZZ, Br_Amumu, Br_Atata, Br_Ass, Br_Abb, Br_Auu, Br_Acc, Br_Agg, Br_Agaga, Br_AgaZ, Br_AZh, Br_AZH, Br_AWHp,
				 gamma_total_Hp, gamma_Hptb, gamma_Hpcs, gamma_Hpcb, gamma_Hptanu, gamma_Hpmunu, gamma_Hpenu, gamma_HpWh, gamma_HpWH, gamma_HpWA,
				                 Br_Hptb, Br_Hpcs, Br_Hpcb, Br_Hptanu, Br_Hpmunu, Br_Hpenu, Br_HpWh, Br_HpWH, Br_HpWA,
				 gamma_total_top, gamma_tHpb, Br_tHpb,
		         CS_singleH, CS_Hhh, CS_bbgaga, CS_tHpb);
	//f.close();
	fclose(fp);
}
